A toy model of polymer stretching 
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We present an extremely simplified model of multiple-domains polymer stretching in an atomic 
force microscopy experiment. We portray each module as a binary set of contacts and decompose 
the system energy into a harmonic term (the cantilever) and long-range interactions terms inside 
each domain. Exact equilibrium computations and Monte Carlo simulations qualitatively reproduce 
the experimental saw-tooth pattern of force-extension profiles, corresponding (in our model) to first- 
order phase transitions. We study the influence of the coupling induced by the cantilever and the 
pulling speed on the relative heights of the force peaks. The results suggest that the increasing 
height of the critical force for subsequent unfolding events is an out-of-equilibrium effect due to a 
finite pulling speed. The dependence of the average unfolding force on the pulling speed is shown 
to reproduce the experimental logarithmic law. 
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I. INTRODUCTION 

The last decade has witnessed a significant advance- 
ment of single molecule manipulation and visualization 
techniques fl, 0, 0, El IS El providing access to the 
distribution of physical properties across many individ- 
ual molecules and not just average properties as was the 
case of traditional biochemical techniques. 

Optical Tweezers (OT) and Atomic Force Microscopy 
(AFM) in particular, enable the study of mechanical 
properties of proteins such as those in the extracellular 
matrix, in the cytoskeleton and in the muscle, that in 
vivo are exposed to stretching forces. The mechanical 
properties of macromolecules obtained in these experi- 
ments may be directly connected to corresponding ther- 
modynamical quantities @ with a bit of caution, since 
mechanical experiments are usually out-of-equilibrium, 
as discussed in this paper. On the other hand, the me- 
chanical properties of biomolecules may be of direct im- 
portance for what concerns bendability (DNA), translo- 
cation of nucleic acids and proteins across cellular mem- 
branes, rigidity and elasticity (structural proteins). 

In a force- measuring AFM experiment, the tip of a mi- 
croscopic cantilever is pressed against a flat, gold-covered 
substrate, coated with a thin layer of protein. Upon re- 
traction of the substrate, a protein molecule that may 
have been adsorbed to the cantilever tip is then stretched. 
Finally, the force the protein opposes to the stretching 
is computed from the cantilever deflection, and a force- 
extension plot is drawn. In a typical stretching exper- 
iment the force-extension curve shows a saw-tooth pat- 
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tern, each peak corresponding to the unfolding of a single 
domain. However, if the protein is composed by several 
different modules, it is difficult to associate each peak 
to the corresponding domain. Moreover, in order to ne- 
glect the tip-substrate interaction, it is necessary to have 
a sufficiently long protein. A common strategy proposed 
in the literature to address these problems is to use an en- 
gineered protein composed by several tandem repeats of 
the same kind of domain. A typical choice is to use do- 
mains belonging to the immunoglobulin fibronectin 
III |irjfl or cadherin superfamilies, characterized by 
/3-sandwich structures. Alpha-helical domains, such as 
those of the cytoskeletal protein spectrin 0], have also 
been used. 

Many data about the mechanical properties of mod- 
ular proteins can be extracted from the force-extension 
profiles. The amplitude of each peak, in fact, represents 
the mechanical stability of the corresponding domain, 
while the spacing between peaks reflects the length of 
the unfolded domain and is thus proportional to the num- 
ber of amino-acids it comprises. The experimental data 
also show a logarithmic dependence of the height of the 
peaks on the pulling velocity, so that higher forces are 
required for domain unraveling when the pulling occurs 
very quickly. 

The features of the force-extension curves (Figure 0) 
are accounted for by the entropic elasticity of polymer 
chains in solution. The entropy of a polymer, in fact, is 
maximal in the random coil state, whereas it tends to 
zero in the fully extended conformation. The force re- 
quired to stretch a polymer thus reflects the entropy loss 
and it grows in a non-linear way as the molecule is length- 
ened. The entropic elasticity of biopolymers is currently 
modeled through the worm-like chain (WLC) |l3| and 
freely-jointed chain (FJC) 01 models. The WLC model 
describes a molecular chain as a deformable rod whose 
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FIG. 1: The force-extension profile produced in a typical 
AFM stretching experiment performed on a modular polypro- 
tein composed by 8 tandem repetitions of a Ig-like domain 
from titin. The construct is terminated by 2 cystein residues 
expressly introduced to form covalent bonds with a gold sur- 
face. The experiment was performed at a constant speed of 
200nm/s. Notice that only 7 unfolding peaks appear in the 
plot because the cantilever tip, by chance, established a con- 
tact with the second domain of the molecule, skipping the 
first one. 



stiffness is determined by the persistence length p (the 
length scale over which the polymer loses orientation or- 
der). The functional relation between the external force 
and the fractional extension z/L [z is the end-to-end dis- 
tance and L is the contour length) is approximatly given 
by the interpolating formula 



k B T 
P 



1 



A{l-z/Lf 



1 

4 



(1) 



where kg is the Boltzmann constant and T is the tem- 
perature. 

In the FJC model, the polymer is portrayed as a chain 
of rigid segments linked by frictionless joints so that each 
segment can point in any direction irrespective of the 
orientations of the others. A measure of the stiffness 
of the chain is represented by the Kuhn segment length 
b (the average length of the segments). The analytic 
relation between the average end-to-end distance (z) and 
the stretching force F is 



(z) = Li coth 
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In our model, the polymer is described by an array of 
binary variables representing native contacts that can be 
in either of two states: formed or broken. The cantilever, 
on the other hand, is modeled as a harmonic spring in 
series with the molecule. The energy of the system is 
the sum of a harmonic term and a term of long-range 
interaction modeling in a bulk, coarse-grained way, the 
chemical interactions stabilizing the native conformation 
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FIG. 2: Force versus extension for various models: WLC (in- 
terpolated formula) with p — 0.4 nm, L = 29 nm, T = 310 K; 
FJC with b = 2p, ISING with a = 1.3p N = L, HOOK (har- 
monic potential) with K — 1.3p. 



of the protein. In fact, rather than providing a detailed 
description, we account for chemical interactions through 
a folding prize attributed to a domain when the fraction 
of intact contacts is above a threshold. We assume any 
contact breakdown to produce an equal increment in the 
molecule length. 

Let us consider first the force-extension characteristic 
of our model in the absence of folding prize. Let F be a 
constant pulling force, n the number of broken contacts 
and a the length increment per cleaved contact. If no 
folding prize is attributed to the molecule in a native-like 
state, the Hamiltonian of the system is H = —Fna and 
the partition function is 
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Notice that this is also the partition function of an Ising- 
like model in one dimension without coupling among 
spins. The average end-to-end distance can thus be com- 
puted as 



(z) = a(n) — 
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The variation of length (Az) versus F in the WLC, 
FJC and our (ISING) models is shown in Figure 
where we have adjusted the constant a (corresponding 
to persistence length p so to make the curve coincide for 
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FIG. 3: Isometric versus isotensional elongations for M = 3, 
A = 1, K = 0.1, N = 10, = 0.2 and beta = 2 (see Section|IT| 
for the illustration of the model). 



small elongations. Despite the extreme simplicity of our 
assumptions, the three curves are qualitatively similar. 
This similarities could be increased by adding contact- 
contact interactions or considering different elongations 
for the various contact breaking events. 

Although the WLC and FJC models more accurately 
represent the physics of a polymer, their statistical me- 
chanics treatment is so complex that one generally em- 
ployes the average formulas (|TJ and @ for each do- 
main, complemented with an independence hypothesis 
of domains and an ad-hoc treatment of the unfolding 
even t ITEl ITrl ITtI ] , based on an extension of Bell's expres- 
sion |18( or full Kramer's theory [19| for the rupture rate 
coefficient in the presence of a time-dependent external 
force. 

The independence assumption is questionable, since all 
domains are coupled by the presence of the cantilever. 
This difference may be explicited by comparying com- 
putations |2fJ in which the position of the cantilever is 
observed (isometric) while the force may fluctuate, with 
computations in which the force is maintained constant 
(isotensional) . We can obtain exact comparisons of the 
two different set-up for our model (see Section [nj , as 
shown in Figure |21 It can be noticed that the isoten- 
sional model does not show any peaks (the peaks are here 
smoothed due to the small length of the single module) . 

Moreover, the peaks in the force-extension profile 
(Fig. ^) are a signature of a first-order transition. As 
we shall show in the following, this transition naturally 
arise in our model due to the long-range coupling (the 
folding prize). 

Despite its extreme simplicity, our approach captures 
some important aspects of the physics of Ig domain 
stretching. Steered molecular dynamics simulations per- 
formed in Schulten's group null, in fact, showed 
that the mechanical unfolding of the 127 module occurs 
only after the breakdown of a patch of six hydrogen bonds 
bridging the A' and G /3-strands. The rupture of these 



critical bonds was shown to be the key event allowing the 
full unraveling of the molecule under an external force. 

Even if this pattern was originally observed in a specific 
protein, it could be hypothesized a more widespread dis- 
tribution. Makarov and coworkers |24j performed Monte 
Carlo simulations of titin forced unfolding. During these 
simulations the number of hydrogen bonds at time t, n(t), 
undergoes a random walk. It was concluded that a crit- 
ical value n# of the number of hydrogen bond does ex- 
ist, such that when n(t) < n&, the domain unfolds very 
rapidly. Makarov also showed that the number of bonds 
is roughly one, when the force is low enough, whereas for 
very large pulling rates (and thus large pulling forces), it 
is likely to be equal to six, recovering the findings by Lu 
and Schulten |23| . 

The accuracy of the model developed by Makarov and 
coworkers, allows quantitative comparisons with experi- 
mental data at the expense of very long simulation times 
and the need to assume the knowledge of the free en- 
ergy profile. Their model is also based on the hypothesis 
of independence of domains which, however, might be 
incompatible with the coupling introduced by the can- 
tilever. Since the transitions shown in AFM experiments 
are out of equilibrium |l5l l25j , thermal fluctuations may 
play a fundamental role. Our model, conversely, is so 
simplified that we can compute exactly the free energy 
of a multi-domain protein for the equilibrium case, and 
perform long simulations in the out of equilibrium case. 

The model we propose shows interesting similarities 
to a Go model with force rescaling. Recent theoretical 
studies ptl |27| show that the ability of Go models to 
simulate the cooperativity of the folding process can be 
enhanced by imparting an extra energetic stabilization to 
the native state so as to simulate specific interactions ap- 
pearing only after the assembly of native-like structures. 
A rigorous approach to simulate the stabilizing interac- 
tions peculiar of the native state, would be the use of 
two different analytic expressions of the force-field inside 
and outside the native basin. The same purpose can be 
pursued through a much simpler strategy 26] , by rescal- 
ing the conformational force when the fraction of native 
contacts crosses a pre-chosen threshold. 

This approach appears to be similar to the one em- 
ployed in our model. The energy function of our model 
comprises an elastic term and a contact term. The har- 
monic term accounts for the elasticity of the cantilever, 
while the protein can be thought of as a soft spring so 
that its contribution to the elasticity of the system is 
negligible. 

The contact term of our energy function, on the other 
hand, is a stabilizing contribution that the protein re- 
ceives only when the fraction of intact contacts crosses 
a threshold. This approach is thus equivalent to a force 
rescaling occurring whenever the polymer enters the na- 
tive basin, and the contact term of the energy function 
appears to be a square well. 

Finally, from a purely formal point of view, our de- 
scription of the polymer as an array of binary variables, 
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FIG. 4: Schematic description of the AFM experimental set- 
up. The polymer is composed by several tandem repeats of 
the same domain in series with a harmonic spring (the can- 
tilever) . 



bears some similarity with the model proposed by Galz- 
itskaya and Finkelstein (GF) |2^. In the GF model, in 
fact, each residue of the polypeptide chain can be either 
in an ordered (native) or disordered (non-native) state, 
encoded by the two possible values of a binary variable. 
This approach, similarly to ours, significantly narrows 
the conformational space, that consists of 2 N conforma- 
tions only, for a polymer with N residues. This approach, 
while drastically reducing the computation time, is in 
agreement with the Zimm-Bragg model |29j . widely em- 
ployed to describe the helix-coil transition in heteropoly- 
mers. 

The outline of the paper is as follows. In Section ITT1 
wc describe the model and the simulation procefdurc; 
in Section IIII Al we study the equilbrium behavio r of a 
single domain; in Section IIII CI we investigate the role 
of fluctuations in the reciprocal influence between suc- 
cessive unfolding events; in Section llll Dl we investigate 
the dependence of the unfolding force on the pulling rate 
and we sketch a simple analytic treatment of the unfold- 
ing probability in the limit of an extremely high pulling 
rate; finally, in Section llVl we draw the conclusions of our 
work. 



II. THE MODEL 

The AFM in its simplest arrangement is just a spring 
(the cantilever) whose deflection, and thus the applied 
force, is measured as a function of its position. The 
system, as shown in Figure can therefore be modeled 
as a harmonic spring in series with a protein composed 
by M tandem repeats of the same domain. Each do- 
main j is simply portrayed as a sequence of contacts 
(sCi) = {sp'}, i = 1, ••• ,JV) that can be either intact 

(s^ = 0) or broken (s^ = 1), where N is the to- 
tal number of native contacts inside each domain. The 
length of the polymer chain can be simply computed as 

I = Yl^Li 12iLi s^ 11 = Sjf=i an ji where a represents the 
incremental elongation associated to each contact break- 
down and rij — ^ s ^ ne number of broken con- 



tacts in domain j. For the sake of simplicity, in all our 
computations we set a = 1. The length of the spring, on 
the other hand, is just x = A — I where A is the extension 
of the spring-polymer system, using the rest position as 
the reference point. 

The energy of a domain configuration s^> is modeled 
by the sum of two contributions: a harmonic term H, 
accounting for the presence of the spring, and the sum 
over all domains j of a term related to the long- 

range interactions among the monomers, 

M 

E = H + J2 lU) ' 

3=1 

The harmonic term is expressed as 

1 1 1 / M V 

H(x) = -Kx 2 = -K(X - if = -K I A - «X> J ' 

where K is the harmonic constant. 

In our simplified representation, if the fraction of intact 
contacts rij/N in a domain j is below a given threshold 

the domain receives a folding prize AN proportional 
to the number of possible contacts. The interaction term 
is thus computed as 

L {j) = -ANQ(9N — rij), 

where Q(x) is the Heavisidc step function 



0(x) 



if x < 

1 otherwise 



In summary, the energy of a configuration is just a 
function of the number of broken contacts in each domain 

E{n) = -K[x(n)f - ANG(9N - rij), 

where n = {n\, n%, . . . , Um} and 
/ 

x{n) = I A — a rij 

This assumption speeds up the computations, that may 
be performed in terms of the rij. 

A stretching simulation starts from a completely folded 
initial structure, where no contact is broken. The 
protein-spring length A, chosen as the control parame- 
ter, is linearly increased from Xmw to \mo,x in k\ + 1 
steps during the simulation; 

\{k) = X M in H : k, 

k\ 

where k = 0, 1, • • • , k\. For each value of A, we compute 
the average length of the spring (a;) in an equilibrium 
simulation as 



(x) = ^J2d( n ) x{n)e 



■pE{n) 
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where Z is the partition function 



and the multeplicity factor 



M 

9(n) = II 

3=1 



is given in terms of the number of possible microscopic 
configurations containing rij cleaved contacts. 

In real stretching experiments, however, the polymer is 
subjected to a finite pulling velocity, so that the molecule 
cannot be considered in an equilibrium condition. In or- 
der to consider this effect, we use Monte Carlo simu- 
lations. For each value of A, T Monte Carlo steps are 
performed, each involving N x M elementary steps. The 
elementary step consists in the random choice of a do- 
main and in the attempt to increase or decrease by one 
the number of contacts with probability equal to the frac- 
tion of broken or intact contacts respectively. The trial 
move is then accepted or rejected with a probability de- 
rived from the heat-bath criterion, 



p(nj -► n 3 ±l) = 



-/3J?(ni,...,nj±l,....nM) 



-0E(ni 
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where (3 is the inverse temperature. 

The average length of the spring (x), corresponding 
to the current position A is computed averaging over T 
Monte Carlo steps. 



III. RESULTS 

We first analyze the entropic effects related to tem- 
perature through the analysis of computations in the ab- 
sence of a folding prize, and then investigate the role 
of long-range interaction by setting a non-zero prize on 
a single-domain polymer. After that, we discuss a set 
of computations on a three-domain protein, showing the 
importance of the coupling due to the cantilever in the 
mechanism of mechanical unfolding and, in particular, 
they explain how the first unfolding event affects the fol- 
lowing ones. In the last part, we discuss the relation 
between pulling rate and unfolding force, finding a log- 
arithmic law. The section is completed with a analytic 
treatment of the unfolding probability valid in the limit 
of high pulling rate. 

The legend of the symbols appearing in the figures is 
shown in Table [i] 



N : 


Maximum number of contacts per domain 


M : 


Number of domains 


L : 


Maximal length of the polymer- spring system 


kx : 


Number of steps in the length of the polymer-spring system 


T : 


Number of Monte Carlo steps 


A : 


Folding prize (in units TV) 


6 : 


Threshold to keep the folding prize 


K : 


Elastic constant of the cantilever spring 


P ■ 


Inverse temperature 



TABLE I: Legend of the symbols appearing in the figures. 




FIG. 5: Cantilever deflection (force) versus estension for a 
single domain, (a) No folding price. Influence of temperature- 
dependent entropic effects on mechanical unfolding. Simula- 
tion parameters: N = 30; M = 1; L = 50; k x = 50; A = 0; 
.= 0.1; P = 0.001 (solid line); p = 50 (dashed line), (b) 



2 typical ramp-like profile with folding prize. Simulation pa- 
rameters: N = 100; L = 300; fc A = 100; A = 5; 6 = 0.5; 
K = 0.1; p = 2. 



part of the curve at low temperature {[3 = 50) represents 
the complete unfolding of the protein while the spring 
nearly retains its resting length: at low temperatures, the 
enthalpic contribution of the free energy (the harmonic 
energy of the cantilever), dominates over the entropic 
one. When the protein is completely stretched, the sys- 
tem can react to the increase of the control parameter A, 
only through an equal increase of the spring length. The 
steep part of the (x, A) plot is thus a straight line with 
unitary slope. 

At high temperature ((3 — 0.001) the free energy is 
dominated by the entropic term so that for small values 
of A, about 50% of the monomers are extended in the 
direction of the pulling force so as to maximize entropy, 
while the spring remains in its resting position. The pro- 
portion of unfolded monomers remains thereafter almost 
unchanged during the simulation and for A > 15, any 
further increase in A is reflected in a equal extension of 
the spring Ax — A A. 



B. Effect of folding prize 



A. Single domain analysis 

The force-extension curves without folding prize, as 
shown in Figure (a), are clearly bi-phasic. The flat 



The protein is here composed by a single domain with 
N = 100 contacts, and its energy is lowered by A = 5N 
when a fraction of residues greater than 9 = 0.5 is folded. 

The (x,A) plot portrayed in Figure [3(b) shows a flat 
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FIG. 6: Role of folding prize A (a) and harmonic constant 
K (b). Common computation parameters: N = 100; M = 1; 
L = 300; k x = 100; 6 = 0.20; /3 = 2. 



region for A < 50. This reflects the cleavage of N x 9 = 50 
contacts that can occur without the loss of the folding 
prize, while the spring remains very close to the resting 
length. As the further extension of the protein would 
result in a significant destabilization of the system due 
to the loss of the folding prize, the increase of the con- 
trol parameter A is now completely accounted for by the 
stretching of the spring. The second part of the (x, A) 
plot is thus a straight line with unit slope. When the in- 
crease in harmonic energy exceeds the folding prize, the 
stretching of the spring is interrupted and the remaining 
50 contacts of the protein break down, allowing a cor- 
responding shortening of the spring. As the protein is 
now completely extended, any further increase in A must 
result in a corresponding stretching of the cantilever and 
the final part of the (x, A) plot is again a straight line 
with unit slope. 

The features of the saw-tooth (x, A) profile are affected 
by several simulation parameters. The folding prize A is 
related to the enthalpic component of the free energy 
and stabilizes the folded conformation of the protein. As 
a consequence, when A is large, the polymer tends to 
remain in the compact conformation so as to retain the 
significant folding prize and the increase in A leads to a 
stretching of the cantilever spring. Only when (x) is very 
large it becomes enthalpically favourable for the polymer 
to unfold because the decrease in harmonic energy due 
to the cantilever relaxation more than compensates the 
loss of the folding prize. Thus, as A is increased, higher 
and higher values of (x) are required for the elastic en- 
ergy to compensate the folding prize and the peak of the 
(x, A) plot becomes accordingly higher and higher. The 
harmonic constant K of the cantilever spring plays a role 
basically opposite to that of the folding prize. In fact, 
when K is large, smaller average extensions (x) are re- 
quired for the harmonic energy to balance the folding 
prize so that larger K s result in a less pronounced peak 
in the (x, A) plot. The role of A and K is exemplified by 
the simulations portrayed in Figure 

As just discussed, the role of the folding prize and of 
the harmonic constant is related to the enthalpic term 
of the free energy. By contrast, the folding threshold 9 
influences the entropic contribution to the free energy. 
A small value of 9 in fact, implies that only a small 
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FIG. 7: Role of the threshold 6. Other computation parame- 
ters: N = 100; M = 1; L = 300; k x = 100; A = 5; K = 0.1; 
13 = 2. 



fraction of the contacts can be broken without loss of 
the folding prize. If m < 9N represents the number 
of broken contacts in a moment preceding the unfold- 
ing event, then the number of microscopic conformations 
allowed will be g{n,N) — (^), and the entropy will be 
S = -T\og[g{n, N)]. If 9 < 0.5, the unraveling of the do- 
main will increase the degeneracy g(n, N) and thus the 
entropy, so that a moderate extension of the spring will 
be sufficient for the enthalpic term to be more than com- 
pensated by the entropic one. By contrast, when 9 > 0.5, 
the breakdown of the domain and the resulting increase 
in the number of disrupted contacts, will bring the degen- 
eracy and the entropy further away from the maximum 
thus disfavouring the unfolding event and causing the 
peak of the (x, A) profile to become higher. The pattern 
of increase in height of the unfolding peak as 9 takes on 
higher values is shown in Figure [3 

The relevance of the entropic contribution on free en- 
ergy computation strongly depends on temperature that 
may amplify the role of the threshold 9. As already no- 
ticed, in fact, when 9 < 0.5, the breakdown of the do- 
main is entropically favoured as it brings the degeneracy 
closer to its maximum. Since this gain in entropy be- 
comes larger and larger as the temperature is increased, 
for small 9 the unfolding event becomes more and more 
favourable as f3 is decreased, and the height of the peak of 
the (x, A) plot will also decrease accordingly. For 9 > 0.5 
the opposite pattern can be observed. In fact, the entropy 
loss due to the decrease of the degeneracy resulting from 
the unfolding event, is magnified as (3 is decreased lead- 
ing to higher and higher peaks in the (x, A) plot. The 
interplay between the parameters 9 and (3 is shown in 
Figure OD 
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FIG. 8: Interplay between the folding threshold 9 and the 
inverse temperature /3. For small 9s (a) higher temperatures 
favour the unfolding, whereas the opposite is true for large 9s 
(b). Common computation parameters: TV = 100; M = 1; 
L = 300; fc A = 100; A = 5; K = 0.1. 
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FIG. 9: Role of fluctuations. Panels (a) and (b): a threshold 
lower than 0.5 (9 = 0.1) causes the first unfolding event to 
favour the following ones. Panels (c) and (d): the first un- 
folding event does not affect the second and third ones as a 
result of the high threshold 8 = 0.5. Computation parame- 
ters: TV = 100; M = 3; L = 400; k x = 400; A = 1; K = 0.05; 
(3 = 2. 



C. Coupling and fluctuations 

Before studying the mutual influence of the unfolding 
events, let us illustrate the features of the free energy 
landscape of a single module near the unfolding transi- 
tion. In Figure El we show the evolution of the free en- 
ergy landscape in correspondence of the unfolding tran- 
sition in a simulation with M = 3 and = 0.1 (see also 
top left panel of Fig.[!||). The profile is characterized by a 
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FIG. 10: Free energy landscape under the action of an in- 
creasing elastic force in correspondence of the unfolding peak 
shown in the top left panel of Figure [5] Simulation parame- 
ters: TV = 100; M = 3; L = 400; fc A = 400; A = 1; 9 = 0.1; 
K = 0.05; [3 = 2. 



cusp-like, narrow well corresponding to the folded state, 
and a wide, smooth well related to the unfolded state. 
The width of the two wells, in fact, depends on the num- 
ber of conformations that the system can explore: in the 
folded state, the system is stretched due to the action of 
the spring and no fluctuations are allowed so that only 
one conformation will be populated; after the transition, 
the residues of the collapsed domain become free to fluc- 
tuate and many conformations will be explored thus de- 
termining a very wide well in the free energy profile. The 
Figure shows that for low values of A (and thus low val- 
ues of the elastic force), the free energy of the reference 
folded state Go is lower than that of the unfolded state 
G, thus forbidding the breakdown of the domain; as A is 
increased, the free energy of the folded state increases, 
until it finally becomes higher than that of the folded 
conformation and the stretching transition occurs. This 
is a typical example of a first-order phase transition. 

We now consider a polymer composed by M = 3 tan- 
dem repeats with TV = 100 contacts each. In Figure |5| 
we compare two simulations performed with the same 
parameter set except for a different threshold 9. In the 
simulation with 9 = 0.1 the peaks corresponding to the 
second and third unfolding events are less pronounced 
than the first one thus suggesting that the unfolding of a 
domain actually favours further unfolding events. Con- 
versely, in the simulation with 9 = 0.5, the three peaks 
feature almost the same height showing that the first 
unfolding event has little or no influence at all on the 
following ones. 

Each unfolding event in Figure [5] corresponds to a 
peak in the variance plot because the unraveling of a do- 
main increases the fluctuations of the polymer and spring 
length. The variance peaks are characterized by a high 
and narrow spike followed by a smoother region that de- 
creases more slowly. The shape of the variance peak is 
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related to the regions of the free energy landscape ex- 
plored by the system during the unfolding transition 
The spikes in the variance plots (that are truncated foi 
the sake of graphical clarity) correspond to the situatior 
with G = Go when both wells are explored by the systen 
and the variance a^. is related to the distance between the 
two wells. On the other hand, the smooth regions of the 
variance plots refer to the case with G < Gq when the 
system only explores the unfolded region of the free en- 
ergy landscape whose width correlates with the variance 

Figure [5] shows that the height of the smooth regior 
of the variance peaks increases with the order of the un- 
folding event. This trend is due to the fact that, with 
each unfolding event, N(l — 9) new monomers become 
free to fluctuate and the number of accessible conforma- 
tions increases accordingly. Figure also shows that the 
value of the variance o\ after each unfolding event, grad- 
ually decrease as A is increased, because the disruption of 
the contacts of the domain just collapsed allow an exten- 
sion of the molecule in the direction of the pulling force 
so as to avoid as far as possible a further stretching of 
the spring that would cause an increase in energy. As 
a result, narrower and narrower regions of the confor- 
mation space become accessible to the polymer and the 
variance is lowered. It is worthwhile noticing, however, 
that the extension of the unfolded domain is hindered by 
the subsequent decrease in entropy so that the number 
of contacts actually broken before each unfolding event 
is smaller than the maximum number allowed by the loss 
of the folding prize in the previous domain breakdown. 

When 9 = 0.1, 26 contacts (a number of the order of 
9NM) are broken before the first unfolding event. This 
value is consistent with the number of contacts that can 
be disrupted without loss of the folding prize. After the 
first collapse event, the number of contacts broken in the 
simulation rises to 115, thus determining an increase of 
the fluctuations and favouring the breakdown of another 
domain. This explains why the second peak of the (x, A) 
plot is less pronounced than the first one. The occurrence 
of the second unfolding determines a further increase of 
the number of cleaved contacts to 205. This results in 
a easy breakdown of the last domain of the protein and 
the last peak of the (x, A) plot is therefore less high than 
the second one. 

This scenario is significantly different for 9 = 0.5. For 
large values of 9, in fact, only a small number of residues 
can be recruited for fluctuations after each unfolding 
event. As a consequence, the variance a^. rapidly goes 
to zero after each unfolding event and the fluctuations 
are extinguished before the force threshold for unfolding 
can be crossed. The following unfolding events, similarly 
to the first one, will occur when the increase in harmonic 
energy balances the loss of the folding prize and therefore 
the height of the unfolding (x)-peaks will be roughly the 
same. 

The role of fluctuations in determining the relative 
heights of the peaks of the (x, A) plot, and thus the cou- 
pling among domains, is confirmed through simulations 
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FIG. 11: Effect of temperature on the relative heights of the 
peaks of the saw-tooth profile. For 9 = 0.1, if the tempera- 
ture is sufficiently high (/3 = 2, panels (a) and (b)), the first 
unfolding event favours the following ones due to the role of 
fluctuations. At low temperatures however (f3 — 8, panels 
(c) and (d) ), the fluctuations rapidly become negligible and 
the heights of the peaks become roughly the same due to the 
identical energetic features of the domains. Computation pa- 
rameters: N = 100; M = 3; L = 400; = 400; A = 1; 
9 = 0.1; K = 0.05. 

performed at different temperatures. At very low tem- 
peratures, in fact, the entropic term of the free energy 
becomes negligible and the height of the peaks of the 
(x, A) plot only depends on the energetic term. As we are 
considering a protein composed by identical domains, we 
expect the x(X) peaks to be identical if the temperature 
is sufficiently low. This pattern can be observed in Fig- 
ure ^2 where we compare two simulations performed at 
different temperatures, namely (3 = 2 and (3 — 8. 

D. Pulling rate effects 

Up to now, we discussed equilibrium stretching com- 
putations, i.e. we assumed an infinitely slow pulling. 
However, at the typical time scales of an AFM stretch- 
ing experiment, the polymer is pulled at a finite velocity 
and the unfolding is a non-equilibrium process, as tes- 
tified by the differences between the unfolding and the 
refolding force-extension profiles (see e.g. Figure IT^ . In 
fact, while the unfolding profile features the typical saw- 
tooth pattern, upon relaxation of the unfolded polymer, 
the trace exhibits no discontinuities that would indicate 
refolding. 

A consequence of the irreversibility of the stretching 
process is that the unfolding force depends on the pulling 
speed. Actually, if the loading rate rj — kfV (with kf be- 
ing the elastic constant and v the pulling velocity) is suf- 
ficiently small, thermal fluctuations are allowed enough 
time to overcome the energy barrier and the unfolding 
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FIG. 12: Monte Carlo computations of stretching and re- 
folding: the non-superposability of the profiles is a signature 
of the irreversibility of the process. Simulation parameters: 
N = 100; M = 3; L = 450; k x = 100; T = 1000; A = 1; 
9 = 0.3; K = 0.5; = 2. 



force will be low. 

Several experimental works 0, H3, El H2 reported 
a logarithmic dependence of the unfolding force on the 
loading rate in the case where a single energy barrier is 
present along the reaction path. The analytic expression 
of the relation between force and loading rate was de- 
rived [Tslls^lls^ l within the frame of Kramer's theory for 
a simple two-state model, by assuming that the external 
force reduces the height of the energy barrier, 



k B T ( KvAx \ 



(4) 



where ks is Boltzmann constant, T is the absolute tem- 
perature, Ax is the distance between the minimum cor- 
responding to the folded state and the activation barrier 
of the energy landscape, v is the pulling speed, K is the 
cantilever harmonic constant and ko is the spontaneous 
unfolding rate. 

Recently 19], it has been shown that the probability 
distribution of the force at various pulling rates does not 
follow the simple Bell law, requiring the full Kramer's 
theory and predicting small corrections to the logaritmic 
behavior of the most probable force. 

We limit here to a preliminary illustration of a series of 
Monte Carlo simulations with different number of steps 
T considering the pulling speed as being proportional to 
l/T. A more detailed analysis will be presented in a 
forthcoming paper. 

Figure FLU shows that, as the number of Monte-Carlo 
(MC) steps is increased, the (x — A) plot becomes closer 
and closer to the profile yielded by the equilibrium simu- 
lation. In particular, a small number of MC steps, i.e. a 
fast pulling, corresponds to high peaks of the (x — X) pro- 
file, whereas larger numbers of steps correspond to less 



FIG. 13: Comparison between the equilibrium computations 
and a set of Monte Carlo simulations with different numbers 
of MC steps. The peaks of the (x — A) profile become lower 
and lower and the MC simulations converge to the equilibrium 
scenario as the number of MC steps is increased. Simulation 
parameters: N = 10; M = 3; L = 70; k x = 100; A = 1; 
8 = 0.2; K = 0.1; = 2. Notice that, using a smaller value of 
N, the equilibrium profile is smoother than that in Figure |5] 
Figure rm and preceding ones. 
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FIG. 14: Statistical analysis of unfolding forces in 5 sets of 
100 independent Monte Carlo runs with different numbers of 
MC steps. Panel (a): histograms of spring extensions at the 
rupture point; panel (b): linear fit of the average rupture force 
as a function of the logarithm of the loading rate. Simulation 
parameters: N = 10; M = 3; L = 70; fc A = 100; A = 10; 
9 = 0.2; K = 0.1; = 2. 



and less pronounced peaks and thus, smaller unfolding 
forces. 

We performed a series of 100 independent MC runs for 
5 different values of T: 1 ■ 10 3 , 5 ■ 10 3 , 1 • 10 4 , 5 • 10 4 and 
1 • 10 5 . For each series of runs we built the histograms of 
the rupture spring elongations and we plotted the mean 
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value of the histogram as a function of the logarithm 
of the loading rate. Finally, the set of points thus ob- 
tained was linearly fitted. As shown in Figured! the his- 
tograms, shift to lower values as the number of MC steps 
increases. Figure 1141 also shows that the mean rupture 
forces computed from the histograms, feature an excel- 
lent linear correlation with the logarithm of the loading 
rate (correlation coefficient r c = 0.99). 

If Eq. 0] is explicitly rearranged so as to show the lin- 
ear dependence of (F max ) on log(Kv), the coefficients of 
the linear equation can be equated to the corresponding 
parameters 71 and 72 of the regression line (F max ) = 
71 log(Kv) +72, so as to build the following system of 
equations: 



k B T 
Ax 



Ax 
k k B T 



71 
72 



From the first equation of the system it is possible to 
compute the width of the activation barrier Ax = 11.17; 
this value can then be substituted into the second equa- 
tion so as to determine k = 2.93 • 10~ 7 . The spontaneous 
unfolding rate fco, on turn, is related to the height of the 
activation barrier for the unfolding process, 



fc 



: toe 



-AG'l/k B T 



where w, as explained by Kramer's theory, is the recipro- 
cal of a diffusive relaxation time. This simple computa- 
tion shows how stretching experiments and simulations 
provide easy access to important features of the free en- 
ergy landscape. In a forthcoming paper we aim at inves- 
tigating the relations relating the barrier width Ax and 
the spontaneous unfolding rate ko with molecular prop- 
erties such as the folding prize A, the threshold 9, the 
number N and the length M of the domains. 



E. High pulling rate limit 

Let us investigate the dependence of the Monte Carlo 
simulations on the number of Monte-Carlo steps T, that 
we can interpret as a measure of the pulling speed. In 
fact, in the limit of an extremely high pulling speed, we 
can assume that for each value of A the polymer can 
adopt just a single (or at most, a few) conformation, so 
that the entropic contribution can be neglected in the 
computations. 

The mechanism outlined in Section IIII Bl shows that 
the key event is the loss of the folding prize determining 
the complete extension of the protein domain and the 
subsequent relaxation of the spring. The problem thus 
arises to identify the factors affecting the probability of 
this crucial step. More specifically, suppose in a domain 
a number of contacts just below the threshold to keep the 
folding prize has been broken. What is the probability tt 
that one more contact will be cleaved within the next v 
steps ? 




FIG. 15: Effect of the simulation parameters on the height of 
the unfolding ramp. Unless otherwise indicated in the legend, 
the simulation parameters are: TV = 100; M = 1; L = 150; 
k x = 50; T = 100; A = 0.1; 9 = 0.2; K = 0.1; (3 = 2. 



In order to answer this question, it is necessary to com- 
pute the energy difference associated with the transition. 
Let x be the length of the spring when a number of con- 
tacts just below the threshold ON has been disrupted: 
the system retains the prize A. If one more contact is 
broken, the prize is lost and the spring will be corre- 
spondingly shortened. In particular, by assuming each 
contact breakdown to cause a unit increase in length of 
the polymer, the new length of the spring will be x — 1. 
The energy difference between the two states can thus be 
computed as 



AE = ( ^Kx 2 



A 



-K(x - l) 2 = Kx - A, 



where a term K/2 has been neglected. We can now com- 
pute the probability to destroy a contact in a single step 



1 



1 



P 



1 



-0AE 



1 + e (i{A=Kx) ■ 



Finally, the probability to break a contact in v steps is 
given by the sum of a geometric progression 

7r=^(i=p)'p=i-(i-pr- 

The computations thus show that the domain unfolds 
with a probability depending on the prize, the cantilever 
stiffness, the temperature and the pulling velocity (that 
in our simulations is related to the number T of Monte 
Carlo steps). 

In order to test our analytical treatment of the unfold- 
ing probability, we performed a series of MC simulations 
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differing from a reference one just for a parameter (Fig- 
ure !15|) . The key parameters characterizing the reference 
simulations are: T = 100, A = 0.1, K = 0.1, = 2. 
In agreement with the analytic computation, the simula- 
tions show that an increase in the number of Monte Carlo 
steps (T — 10 4 ), an increase in temperature (f3 — 1.5), 
a decrease of the folding prize (A = 8.5 x 10~ 2 ), and an 
increase of the spring constant (K — 1.25 x 10 _1 ), all 
result in an increase of the unfolding probability, thus 
reducing the height of the peak in the (x — A)-plot. 

IV. CONCLUSIONS 

We developed an extremely simplified model of poly- 
mer stretching in which the molecule is portrayed as a 
series of modules, represented as an array of contacts, and 
a harmonic spring (the cantilever). The chemical in- 
teractions stabilizing the native conformation are sim- 
ply modeled as a folding prize gained by domains where 
the fraction of folded monomers is above a pre-chosen 
threshold. Our model is consistent with recent findings 
in Refs [2lL I22I I23I ] , showing that the unraveling of the 
titin Immunoglobulin domain occurs very rapidly only 
after the breakdown of a critical number of key hydrogen- 
bonds. The attribution of a folding prize when a thresh- 
old value of the fraction of native contacts is crossed, also 
makes our approach equivalent to a Go-model with force 
rescaling |26j. However, our model is thus significantly 
simpler than other models commonly used to study me- 
chanical unfolding such as the WLC and FJC models 
and the detailed all-atom, topological and united-residue 
models employed for steered molecular dynamics simu- 
lations. Yet, our model is detailed enough to reproduce 
many qualitative features of the force-extension profiles 
recorded in AFM experiments. In particular, our model 
correctly reproduces the typical saw-tooth pattern with 
peaks characterized by a height dependent on the folding 
prize, the temperature and the pulling velocity. 

In our study, a particular attention was paid to the 
relation between the heights of successive peaks in the 
force-extension plot. This point is quite intriguing as ex- 
perimental data show that the force of unfolding tends 
to increase with each unfolding event, suggesting that 
the protein domains unfold following an increasing order 
of mechanical stability j^Ej. A possible explanation sug- 
gested by the literature is that the domains of the engi- 
neered modular constructs used in AFM experiments are 
not identical but just structurally similar 15] . This ar- 
gument may be correct in the case of real proteins where 
the differences in mechanical stability of the tandem re- 
peats of the construct may arise from their different po- 
sition in the tridimensional structure of the molecule. In 



a computer simulation, however, the domains are all per- 
fectly identical and the explanation for the staircase pat- 
tern of unfolding peaks must be sought elsewhere. This 
behavior has been ascribed to an independent breaking 
probability of each domain. This hypothesis is however 
not consistent with the fact that the cantilever couples 
the fluctuations of all domains. Our unified treatment 
allows to keep into consideration all contributions, and 
to obtain the correct equilibrium curves, that generally 
do not exhibit this effect. However, a detailed study of 
out-of-cquilibrium extension curves, shows that this be- 
havior is consistent with a finite pulling speed, i.e. it can 
be ascribed to a dynamical source. 

Our results appear to be in agreement with recent find- 
ings by Cieplak 36] in steered molecular dynamics simu- 
lations of titin and calmodulin unfolding using a Go-like 
model. In particular, it was found that an increase in ther- 
mal fluctuations results in a lowering of the force peaks 
and in their earlier occurrence during stretching. 

An interesting feature of AFM stretching experiments 
is that the mean unfolding force is a logarithmic function 
of the pulling rate. This pattern reflects the role of ther- 
mal fluctuations: if the pulling speed is sufficiently low, 
then there will be enough time for thermal fluctuations 
to drive the polymer over the free energy barrier and the 
unfolding force will be accordingly low. The ability to 
reproduce this pattern constitutes a benchmark for any 
stretching model. In order to test our model, we per- 
formed a series of Monte Carlo simulations with different 
numbers T of MC-steps and regarding the pulling veloc- 
ity to be proportional to 1/T. Our toy-model, despite its 
extreme simplicity, correctly reproduced the linear de- 
pendence of the mean unfolding force on the logarithm 
of the pulling rate. 

In summary, our simple model including harmonic and 
long-range energy contributions, is not only capable of re- 
producing the force-extension saw-tooth pattern, but it 
also yields force spectra displaying the correct logarith- 
mic dependence of (F max ) on Kv reported in experimen- 
tal works. Finally, our model provides a simple expla- 
nation of the influence of each unfolding event on the 
following ones, based on the role of fluctuations. Our 
work thus shows how minimal models can be valuable 
tools even in the study of complex molecular systems. 
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We present an extremely simplified model of multiple-domains polymer stretching in an atomic 
force microscopy experiment. We portray each module as a binary set of contacts and decompose 
the system energy into a harmonic term (the cantilever) and long-range interactions terms inside 
each domain. Exact equilibrium computations and Monte Carlo simulations qualitatively reproduce 
the experimental saw-tooth pattern of force-extension profiles, corresponding (in our model) to first- 
order phase transitions. We study the influence of the coupling induced by the cantilever and the 
pulling speed on the relative heights of the force peaks. The results suggest that the increasing 
height of the critical force for subsequent unfolding events is an out-of-equilibrium effect due to a 
finite pulling speed. The dependence of the average unfolding force on the pulling speed is shown 
to reproduce the experimental logarithmic law. 
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I. INTRODUCTION 

The last decade has witnessed a significant advance- 
ment of single molecule manipulation and visualization 
techniques fl, 0, 0, El IS El providing access to the 
distribution of physical properties across many individ- 
ual molecules and not just average properties as was the 
case of traditional biochemical techniques. 

Optical Tweezers (OT) and Atomic Force Microscopy 
(AFM) in particular, enable the study of mechanical 
properties of proteins such as those in the extracellular 
matrix, in the cytoskeleton and in the muscle, that in 
vivo are exposed to stretching forces. The mechanical 
properties of macromolecules obtained in these experi- 
ments may be directly connected to corresponding ther- 
modynamical quantities @ with a bit of caution, since 
mechanical experiments are usually out-of-equilibrium, 
as discussed in this paper. On the other hand, the me- 
chanical properties of biomolecules may be of direct im- 
portance for what concerns bendability (DNA), translo- 
cation of nucleic acids and proteins across cellular mem- 
branes, rigidity and elasticity (structural proteins). 

In a force- measuring AFM experiment, the tip of a mi- 
croscopic cantilever is pressed against a flat, gold-covered 
substrate, coated with a thin layer of protein. Upon re- 
traction of the substrate, a protein molecule that may 
have been adsorbed to the cantilever tip is then stretched. 
Finally, the force the protein opposes to the stretching 
is computed from the cantilever deflection, and a force- 
extension plot is drawn. In a typical stretching exper- 
iment the force-extension curve shows a saw-tooth pat- 
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tern, each peak corresponding to the unfolding of a single 
domain. However, if the protein is composed by several 
different modules, it is difficult to associate each peak 
to the corresponding domain. Moreover, in order to ne- 
glect the tip-substrate interaction, it is necessary to have 
a sufficiently long protein. A common strategy proposed 
in the literature to address these problems is to use an en- 
gineered protein composed by several tandem repeats of 
the same kind of domain. A typical choice is to use do- 
mains belonging to the immunoglobulin fibronectin 
III |irjfl or cadherin superfamilies, characterized by 
/3-sandwich structures. Alpha-helical domains, such as 
those of the cytoskeletal protein spectrin 0], have also 
been used. 

Many data about the mechanical properties of mod- 
ular proteins can be extracted from the force-extension 
profiles. The amplitude of each peak, in fact, represents 
the mechanical stability of the corresponding domain, 
while the spacing between peaks reflects the length of 
the unfolded domain and is thus proportional to the num- 
ber of amino-acids it comprises. The experimental data 
also show a logarithmic dependence of the height of the 
peaks on the pulling velocity, so that higher forces are 
required for domain unraveling when the pulling occurs 
very quickly. 

The features of the force-extension curves (Figure 0) 
are accounted for by the entropic elasticity of polymer 
chains in solution. The entropy of a polymer, in fact, is 
maximal in the random coil state, whereas it tends to 
zero in the fully extended conformation. The force re- 
quired to stretch a polymer thus reflects the entropy loss 
and it grows in a non-linear way as the molecule is length- 
ened. The entropic elasticity of biopolymers is currently 
modeled through the worm-like chain (WLC) |l3| and 
freely-jointed chain (FJC) 01 models. The WLC model 
describes a molecular chain as a deformable rod whose 
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FIG. 1: The force-extension profile produced in a typical 
AFM stretching experiment performed on a modular polypro- 
tein composed by 8 tandem repetitions of a Ig-like domain 
from titin. The construct is terminated by 2 cystein residues 
expressly introduced to form covalent bonds with a gold sur- 
face. The experiment was performed at a constant speed of 
200nm/s. Notice that only 7 unfolding peaks appear in the 
plot because the cantilever tip, by chance, established a con- 
tact with the second domain of the molecule, skipping the 
first one. 



stiffness is determined by the persistence length p (the 
length scale over which the polymer loses orientation or- 
der). The functional relation between the external force 
and the fractional extension z/L [z is the end-to-end dis- 
tance and L is the contour length) is approximatly given 
by the interpolating formula 
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where kg is the Boltzmann constant and T is the tem- 
perature. 

In the FJC model, the polymer is portrayed as a chain 
of rigid segments linked by frictionless joints so that each 
segment can point in any direction irrespective of the 
orientations of the others. A measure of the stiffness 
of the chain is represented by the Kuhn segment length 
b (the average length of the segments). The analytic 
relation between the average end-to-end distance (z) and 
the stretching force F is 



(z) = Li coth 
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In our model, the polymer is described by an array of 
binary variables representing native contacts that can be 
in either of two states: formed or broken. The cantilever, 
on the other hand, is modeled as a harmonic spring in 
series with the molecule. The energy of the system is 
the sum of a harmonic term and a term of long-range 
interaction modeling in a bulk, coarse-grained way, the 
chemical interactions stabilizing the native conformation 
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FIG. 2: Force versus extension for various models: WLC (in- 
terpolated formula) with p — 0.4 nm, L = 29 nm, T = 310 K; 
FJC with b = 2p, ISING with a = 1.3p N = L, HOOK (har- 
monic potential) with K — 1.3p. 



of the protein. In fact, rather than providing a detailed 
description, we account for chemical interactions through 
a folding prize attributed to a domain when the fraction 
of intact contacts is above a threshold. We assume any 
contact breakdown to produce an equal increment in the 
molecule length. 

Let us consider first the force-extension characteristic 
of our model in the absence of folding prize. Let F be a 
constant pulling force, n the number of broken contacts 
and a the length increment per cleaved contact. If no 
folding prize is attributed to the molecule in a native-like 
state, the Hamiltonian of the system is H = —Fna and 
the partition function is 
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Notice that this is also the partition function of an Ising- 
like model in one dimension without coupling among 
spins. The average end-to-end distance can thus be com- 
puted as 
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The variation of length (Az) versus F in the WLC, 
FJC and our (ISING) models is shown in Figure 
where we have adjusted the constant a (corresponding 
to persistence length p so to make the curve coincide for 
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FIG. 3: Isometric versus isotensional elongations for M = 3, 
A = 1, K = 0.1, N = 10, = 0.2 and beta = 2 (see Section|IT| 
for the illustration of the model). 



small elongations. Despite the extreme simplicity of our 
assumptions, the three curves are qualitatively similar. 
This similarities could be increased by adding contact- 
contact interactions or considering different elongations 
for the various contact breaking events. 

Although the WLC and FJC models more accurately 
represent the physics of a polymer, their statistical me- 
chanics treatment is so complex that one generally em- 
ployes the average formulas (|TJ and @ for each do- 
main, complemented with an independence hypothesis 
of domains and an ad-hoc treatment of the unfolding 
even t ITEl ITrl ITtI ] , based on an extension of Bell's expres- 
sion |18( or full Kramer's theory [19| for the rupture rate 
coefficient in the presence of a time-dependent external 
force. 

The independence assumption is questionable, since all 
domains are coupled by the presence of the cantilever. 
This difference may be explicited by comparying com- 
putations |2fJ in which the position of the cantilever is 
observed (isometric) while the force may fluctuate, with 
computations in which the force is maintained constant 
(isotensional) . We can obtain exact comparisons of the 
two different set-up for our model (see Section [nj , as 
shown in Figure |21 It can be noticed that the isoten- 
sional model does not show any peaks (the peaks are here 
smoothed due to the small length of the single module) . 

Moreover, the peaks in the force-extension profile 
(Fig. ^) are a signature of a first-order transition. As 
we shall show in the following, this transition naturally 
arise in our model due to the long-range coupling (the 
folding prize). 

Despite its extreme simplicity, our approach captures 
some important aspects of the physics of Ig domain 
stretching. Steered molecular dynamics simulations per- 
formed in Schulten's group null, in fact, showed 
that the mechanical unfolding of the 127 module occurs 
only after the breakdown of a patch of six hydrogen bonds 
bridging the A' and G /3-strands. The rupture of these 



critical bonds was shown to be the key event allowing the 
full unraveling of the molecule under an external force. 

Even if this pattern was originally observed in a specific 
protein, it could be hypothesized a more widespread dis- 
tribution. Makarov and coworkers |24j performed Monte 
Carlo simulations of titin forced unfolding. During these 
simulations the number of hydrogen bonds at time t, n(t), 
undergoes a random walk. It was concluded that a crit- 
ical value n# of the number of hydrogen bond does ex- 
ist, such that when n(t) < n&, the domain unfolds very 
rapidly. Makarov also showed that the number of bonds 
is roughly one, when the force is low enough, whereas for 
very large pulling rates (and thus large pulling forces), it 
is likely to be equal to six, recovering the findings by Lu 
and Schulten |23| . 

The accuracy of the model developed by Makarov and 
coworkers, allows quantitative comparisons with experi- 
mental data at the expense of very long simulation times 
and the need to assume the knowledge of the free en- 
ergy profile. Their model is also based on the hypothesis 
of independence of domains which, however, might be 
incompatible with the coupling introduced by the can- 
tilever. Since the transitions shown in AFM experiments 
are out of equilibrium |l5l l25j , thermal fluctuations may 
play a fundamental role. Our model, conversely, is so 
simplified that we can compute exactly the free energy 
of a multi-domain protein for the equilibrium case, and 
perform long simulations in the out of equilibrium case. 

The model we propose shows interesting similarities 
to a Go model with force rescaling. Recent theoretical 
studies ptl |27| show that the ability of Go models to 
simulate the cooperativity of the folding process can be 
enhanced by imparting an extra energetic stabilization to 
the native state so as to simulate specific interactions ap- 
pearing only after the assembly of native-like structures. 
A rigorous approach to simulate the stabilizing interac- 
tions peculiar of the native state, would be the use of 
two different analytic expressions of the force-field inside 
and outside the native basin. The same purpose can be 
pursued through a much simpler strategy 26] , by rescal- 
ing the conformational force when the fraction of native 
contacts crosses a pre-chosen threshold. 

This approach appears to be similar to the one em- 
ployed in our model. The energy function of our model 
comprises an elastic term and a contact term. The har- 
monic term accounts for the elasticity of the cantilever, 
while the protein can be thought of as a soft spring so 
that its contribution to the elasticity of the system is 
negligible. 

The contact term of our energy function, on the other 
hand, is a stabilizing contribution that the protein re- 
ceives only when the fraction of intact contacts crosses 
a threshold. This approach is thus equivalent to a force 
rescaling occurring whenever the polymer enters the na- 
tive basin, and the contact term of the energy function 
appears to be a square well. 

Finally, from a purely formal point of view, our de- 
scription of the polymer as an array of binary variables, 
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bears some similarity with the model proposed by Galz- 
itskaya and Finkelstein (GF) [2g. In the GF model, in 
fact, each residue of the polypeptide chain can be either 
in an ordered (native) or disordered (non-native) state, 
encoded by the two possible values of a binary variable. 
This approach, similarly to ours, significantly narrows 
the conformational space, that consists of 2 N conforma- 
tions only, for a polymer with N residues. This approach, 
while drastically reducing the computation time, is in 
agreement with the Zimm-Bragg model |29l . widely em- 
ployed to describe the helix-coil transition in hetcropoly- 
mers. 

The outline of the paper is as follows. In Section ITT1 
we describe the model and the simulation procefdurc; 
in Section IIII Al we study the equilbrium behavio r of a 
single domain; in Section IIII CI we investigate the role 
of fluctuations in the reciprocal influence between suc- 
cessive unfolding events; in Section IlII Dl we investigate 
the dependence of the unfolding force on the pulling rate 
and we sketch a simple analytic treatment of the unfold- 
ing probability in the limit of an extremely high pulling 
rate; finally, in Section llVl we draw the conclusions of our 
work. 



II. THE MODEL 
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FIG. 4: Schematic description of the AFM experimental set- 
up. The polymer is composed by several tandem repeats of 
the same domain in series with a harmonic spring (the can- 
tilever) . 

The AFM in its simplest arrangement is just a spring 
(the cantilever) whose deflection, and thus the applied 
force, is measured as a function of its position. The 
system, as shown in Figure can therefore be modeled 
as a harmonic spring in series with a protein composed 
by M tandem repeats of the same domain. Each do- 
main j is simply portrayed as a sequence of contacts 
(gO) = {sp'}, i — 1, • • ■ , N) that can be either intact 

(s\ — 0) or broken (s^ = 1), where N is the to- 
tal number of native contacts inside each domain. The 
length of the polymer chain can be simply computed as 



N M. 



^2jLi an j7 where a represents the 
incremental elongation associated to each contact break- 



tacts in domain j. For the sake of simplicity, in all our 
computations we set a = 1. The length of the spring, on 
the other hand, is just x = A — I where A is the extension 
of the spring-polymer system, using the rest position as 
the reference point. 

The energy of a domain configuration s"' is modeled 
by the sum of two contributions: a harmonic term H, 
accounting for the presence of the spring, and the sum 
over all domains j of a term Lw, related to the long- 
range interactions among the monomers, 

M 

E = H + J2 lU) ' 

3=1 

The harmonic term is expressed as 

1 1 1 / M V 

H{x) = -Kx 2 = -K{\-lf = -K [X-ajyiA , 

where K is the harmonic constant. 

In our simplified representation, if the fraction of intact 
contacts rij/N in a domain j is below a given threshold 

the domain receives a folding prize AN proportional 
to the number of possible contacts. The interaction term 
is thus computed as 

L {j) = -ANQ(9N — rij), 

where Q(x) is the Heaviside step function 



0(x) 



if x < 

1 otherwise 



In summary, the energy of a configuration is just a 
function of the number of broken contacts in each domain 

E{n) = -K[x(n)f - ANG(9N - rij), 

where n = {n\, n%, . . . , Um} and 
/ 

x{n) = I A — a rij 

This assumption speeds up the computations, that may 
be performed in terms of the rij. 

A stretching simulation starts from a completely folded 
initial structure, where no contact is broken. The 
protein-spring length A, chosen as the control parame- 
ter, is linearly increased from Xmw to \mo,x in k\ + 1 
steps during the simulation; 

\{k) = X M in H : k, 

k\ 

where k = 0, 1, • • • , k\. For each value of A, we compute 
the average length of the spring (a;) in an equilibrium 
simulation as 
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where Z is the partition function 



and the multeplicity factor 



M 

9(n) = II 

3=1 



is given in terms of the number of possible microscopic 
configurations containing rij cleaved contacts. 

In real stretching experiments, however, the polymer is 
subjected to a finite pulling velocity, so that the molecule 
cannot be considered in an equilibrium condition. In or- 
der to consider this effect, we use Monte Carlo simu- 
lations. For each value of A, T Monte Carlo steps are 
performed, each involving N x M elementary steps. The 
elementary step consists in the random choice of a do- 
main and in the attempt to increase or decrease by one 
the number of contacts with probability equal to the frac- 
tion of broken or intact contacts respectively. The trial 
move is then accepted or rejected with a probability de- 
rived from the heat-bath criterion, 



p(nj -► n 3 ±l) = 
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where (3 is the inverse temperature. 

The average length of the spring (x), corresponding 
to the current position A is computed averaging over T 
Monte Carlo steps. 



III. RESULTS 

We first analyze the entropic effects related to tem- 
perature through the analysis of computations in the ab- 
sence of a folding prize, and then investigate the role 
of long-range interaction by setting a non-zero prize on 
a single-domain polymer. After that, we discuss a set 
of computations on a three-domain protein, showing the 
importance of the coupling due to the cantilever in the 
mechanism of mechanical unfolding and, in particular, 
they explain how the first unfolding event affects the fol- 
lowing ones. In the last part, we discuss the relation 
between pulling rate and unfolding force, finding a log- 
arithmic law. The section is completed with a analytic 
treatment of the unfolding probability valid in the limit 
of high pulling rate. 

The legend of the symbols appearing in the figures is 
shown in Table [i] 



N : 


Maximum number of contacts per domain 


M : 


Number of domains 


L : 


Maximal length of the polymer- spring system 


kx : 


Number of steps in the length of the polymer-spring system 


T : 


Number of Monte Carlo steps 


A : 


Folding prize (in units TV) 


6 : 


Threshold to keep the folding prize 


K : 


Elastic constant of the cantilever spring 


P ■ 


Inverse temperature 



TABLE I: Legend of the symbols appearing in the figures. 




FIG. 5: Cantilever deflection (force) versus estension for a 
single domain, (a) No folding price. Influence of temperature- 
dependent entropic effects on mechanical unfolding. Simula- 
tion parameters: N = 30; M = 1; L = 50; k x = 50; A = 0; 
.= 0.1; P = 0.001 (solid line); p = 50 (dashed line), (b) 



2 typical ramp-like profile with folding prize. Simulation pa- 
rameters: N = 100; L = 300; fc A = 100; A = 5; 6 = 0.5; 
K = 0.1; p = 2. 



part of the curve at low temperature {[3 = 50) represents 
the complete unfolding of the protein while the spring 
nearly retains its resting length: at low temperatures, the 
enthalpic contribution of the free energy (the harmonic 
energy of the cantilever), dominates over the entropic 
one. When the protein is completely stretched, the sys- 
tem can react to the increase of the control parameter A, 
only through an equal increase of the spring length. The 
steep part of the (x, A) plot is thus a straight line with 
unitary slope. 

At high temperature ((3 — 0.001) the free energy is 
dominated by the entropic term so that for small values 
of A, about 50% of the monomers are extended in the 
direction of the pulling force so as to maximize entropy, 
while the spring remains in its resting position. The pro- 
portion of unfolded monomers remains thereafter almost 
unchanged during the simulation and for A > 15, any 
further increase in A is reflected in a equal extension of 
the spring Ax — A A. 



B. Effect of folding prize 



A. Single domain analysis 

The force-extension curves without folding prize, as 
shown in Figure (a), are clearly bi-phasic. The flat 



The protein is here composed by a single domain with 
N = 100 contacts, and its energy is lowered by A = 5N 
when a fraction of residues greater than 9 = 0.5 is folded. 

The (x,A) plot portrayed in Figure [3(b) shows a flat 



200 





100 150 200 250 300 50 100 150 200 250 300 



FIG. 6: Role of folding prize A (a) and harmonic constant 
K (b). Common computation parameters: N = 100; M = 1; 
L = 300; k x = 100; 6 = 0.20; /3 = 2. 



region for A < 50. This reflects the cleavage of N x 9 = 50 
contacts that can occur without the loss of the folding 
prize, while the spring remains very close to the resting 
length. As the further extension of the protein would 
result in a significant destabilization of the system due 
to the loss of the folding prize, the increase of the con- 
trol parameter A is now completely accounted for by the 
stretching of the spring. The second part of the (x, A) 
plot is thus a straight line with unit slope. When the in- 
crease in harmonic energy exceeds the folding prize, the 
stretching of the spring is interrupted and the remaining 
50 contacts of the protein break down, allowing a cor- 
responding shortening of the spring. As the protein is 
now completely extended, any further increase in A must 
result in a corresponding stretching of the cantilever and 
the final part of the (x, A) plot is again a straight line 
with unit slope. 

The features of the saw-tooth (x, A) profile are affected 
by several simulation parameters. The folding prize A is 
related to the enthalpic component of the free energy 
and stabilizes the folded conformation of the protein. As 
a consequence, when A is large, the polymer tends to 
remain in the compact conformation so as to retain the 
significant folding prize and the increase in A leads to a 
stretching of the cantilever spring. Only when (x) is very 
large it becomes enthalpically favourable for the polymer 
to unfold because the decrease in harmonic energy due 
to the cantilever relaxation more than compensates the 
loss of the folding prize. Thus, as A is increased, higher 
and higher values of (x) are required for the elastic en- 
ergy to compensate the folding prize and the peak of the 
(x, A) plot becomes accordingly higher and higher. The 
harmonic constant K of the cantilever spring plays a role 
basically opposite to that of the folding prize. In fact, 
when K is large, smaller average extensions (x) are re- 
quired for the harmonic energy to balance the folding 
prize so that larger K s result in a less pronounced peak 
in the (x, A) plot. The role of A and K is exemplified by 
the simulations portrayed in Figure 

As just discussed, the role of the folding prize and of 
the harmonic constant is related to the enthalpic term 
of the free energy. By contrast, the folding threshold 9 
influences the entropic contribution to the free energy. 
A small value of 9 in fact, implies that only a small 
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FIG. 7: Role of the threshold 6. Other computation parame- 
ters: N = 100; M = 1; L = 300; k x = 100; A = 5; K = 0.1; 
13 = 2. 



fraction of the contacts can be broken without loss of 
the folding prize. If m < 9N represents the number 
of broken contacts in a moment preceding the unfold- 
ing event, then the number of microscopic conformations 
allowed will be g{n,N) — (^), and the entropy will be 
S = -T\og[g{n, N)]. If 9 < 0.5, the unraveling of the do- 
main will increase the degeneracy g(n, N) and thus the 
entropy, so that a moderate extension of the spring will 
be sufficient for the enthalpic term to be more than com- 
pensated by the entropic one. By contrast, when 9 > 0.5, 
the breakdown of the domain and the resulting increase 
in the number of disrupted contacts, will bring the degen- 
eracy and the entropy further away from the maximum 
thus disfavouring the unfolding event and causing the 
peak of the (x, A) profile to become higher. The pattern 
of increase in height of the unfolding peak as 9 takes on 
higher values is shown in Figure [3 

The relevance of the entropic contribution on free en- 
ergy computation strongly depends on temperature that 
may amplify the role of the threshold 9. As already no- 
ticed, in fact, when 9 < 0.5, the breakdown of the do- 
main is entropically favoured as it brings the degeneracy 
closer to its maximum. Since this gain in entropy be- 
comes larger and larger as the temperature is increased, 
for small 9 the unfolding event becomes more and more 
favourable as f3 is decreased, and the height of the peak of 
the (x, A) plot will also decrease accordingly. For 9 > 0.5 
the opposite pattern can be observed. In fact, the entropy 
loss due to the decrease of the degeneracy resulting from 
the unfolding event, is magnified as (3 is decreased lead- 
ing to higher and higher peaks in the (x, A) plot. The 
interplay between the parameters 9 and (3 is shown in 
Figure OD 
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FIG. 8: Interplay between the folding threshold 9 and the 
inverse temperature /3. For small 9s (a) higher temperatures 
favour the unfolding, whereas the opposite is true for large 9s 
(b). Common computation parameters: TV = 100; M = 1; 
L = 300; fc A = 100; A = 5; K = 0.1. 
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FIG. 9: Role of fluctuations. Panels (a) and (b): a threshold 
lower than 0.5 (9 = 0.1) causes the first unfolding event to 
favour the following ones. Panels (c) and (d): the first un- 
folding event does not affect the second and third ones as a 
result of the high threshold 8 = 0.5. Computation parame- 
ters: TV = 100; M = 3; L = 400; k x = 400; A = 1; K = 0.05; 
(3 = 2. 



C. Coupling and fluctuations 

Before studying the mutual influence of the unfolding 
events, let us illustrate the features of the free energy 
landscape of a single module near the unfolding transi- 
tion. In Figure El we show the evolution of the free en- 
ergy landscape in correspondence of the unfolding tran- 
sition in a simulation with M = 3 and = 0.1 (see also 
top left panel of Fig.[!||). The profile is characterized by a 
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FIG. 10: Free energy landscape under the action of an in- 
creasing elastic force in correspondence of the unfolding peak 
shown in the top left panel of Figure [5] Simulation parame- 
ters: TV = 100; M = 3; L = 400; fc A = 400; A = 1; 9 = 0.1; 
K = 0.05; [3 = 2. 



cusp-like, narrow well corresponding to the folded state, 
and a wide, smooth well related to the unfolded state. 
The width of the two wells, in fact, depends on the num- 
ber of conformations that the system can explore: in the 
folded state, the system is stretched due to the action of 
the spring and no fluctuations are allowed so that only 
one conformation will be populated; after the transition, 
the residues of the collapsed domain become free to fluc- 
tuate and many conformations will be explored thus de- 
termining a very wide well in the free energy profile. The 
Figure shows that for low values of A (and thus low val- 
ues of the elastic force), the free energy of the reference 
folded state Go is lower than that of the unfolded state 
G, thus forbidding the breakdown of the domain; as A is 
increased, the free energy of the folded state increases, 
until it finally becomes higher than that of the folded 
conformation and the stretching transition occurs. This 
is a typical example of a first-order phase transition. 

We now consider a polymer composed by M = 3 tan- 
dem repeats with TV = 100 contacts each. In Figure |5| 
we compare two simulations performed with the same 
parameter set except for a different threshold 9. In the 
simulation with 9 = 0.1 the peaks corresponding to the 
second and third unfolding events are less pronounced 
than the first one thus suggesting that the unfolding of a 
domain actually favours further unfolding events. Con- 
versely, in the simulation with 9 = 0.5, the three peaks 
feature almost the same height showing that the first 
unfolding event has little or no influence at all on the 
following ones. 

Each unfolding event in Figure [5] corresponds to a 
peak in the variance plot because the unraveling of a do- 
main increases the fluctuations of the polymer and spring 
length. The variance peaks are characterized by a high 
and narrow spike followed by a smoother region that de- 
creases more slowly. The shape of the variance peak is 
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related to the regions of the free energy landscape ex- 
plored by the system during the unfolding transition 
The spikes in the variance plots (that are truncated foi 
the sake of graphical clarity) correspond to the situatior 
with G = Go when both wells are explored by the systen 
and the variance a^. is related to the distance between the 
two wells. On the other hand, the smooth regions of the 
variance plots refer to the case with G < Gq when the 
system only explores the unfolded region of the free en- 
ergy landscape whose width correlates with the variance 

Figure [5] shows that the height of the smooth regior 
of the variance peaks increases with the order of the un- 
folding event. This trend is due to the fact that, with 
each unfolding event, N(l — 9) new monomers become 
free to fluctuate and the number of accessible conforma- 
tions increases accordingly. Figure also shows that the 
value of the variance o\ after each unfolding event, grad- 
ually decrease as A is increased, because the disruption of 
the contacts of the domain just collapsed allow an exten- 
sion of the molecule in the direction of the pulling force 
so as to avoid as far as possible a further stretching of 
the spring that would cause an increase in energy. As 
a result, narrower and narrower regions of the confor- 
mation space become accessible to the polymer and the 
variance is lowered. It is worthwhile noticing, however, 
that the extension of the unfolded domain is hindered by 
the subsequent decrease in entropy so that the number 
of contacts actually broken before each unfolding event 
is smaller than the maximum number allowed by the loss 
of the folding prize in the previous domain breakdown. 

When 9 = 0.1, 26 contacts (a number of the order of 
9NM) are broken before the first unfolding event. This 
value is consistent with the number of contacts that can 
be disrupted without loss of the folding prize. After the 
first collapse event, the number of contacts broken in the 
simulation rises to 115, thus determining an increase of 
the fluctuations and favouring the breakdown of another 
domain. This explains why the second peak of the (x, A) 
plot is less pronounced than the first one. The occurrence 
of the second unfolding determines a further increase of 
the number of cleaved contacts to 205. This results in 
a easy breakdown of the last domain of the protein and 
the last peak of the (x, A) plot is therefore less high than 
the second one. 

This scenario is significantly different for 9 = 0.5. For 
large values of 9, in fact, only a small number of residues 
can be recruited for fluctuations after each unfolding 
event. As a consequence, the variance a^. rapidly goes 
to zero after each unfolding event and the fluctuations 
are extinguished before the force threshold for unfolding 
can be crossed. The following unfolding events, similarly 
to the first one, will occur when the increase in harmonic 
energy balances the loss of the folding prize and therefore 
the height of the unfolding (x)-peaks will be roughly the 
same. 

The role of fluctuations in determining the relative 
heights of the peaks of the (x, A) plot, and thus the cou- 
pling among domains, is confirmed through simulations 
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FIG. 11: Effect of temperature on the relative heights of the 
peaks of the saw-tooth profile. For 9 = 0.1, if the tempera- 
ture is sufficiently high (/3 = 2, panels (a) and (b)), the first 
unfolding event favours the following ones due to the role of 
fluctuations. At low temperatures however (f3 — 8, panels 
(c) and (d) ), the fluctuations rapidly become negligible and 
the heights of the peaks become roughly the same due to the 
identical energetic features of the domains. Computation pa- 
rameters: N = 100; M = 3; L = 400; = 400; A = 1; 
9 = 0.1; K = 0.05. 

performed at different temperatures. At very low tem- 
peratures, in fact, the entropic term of the free energy 
becomes negligible and the height of the peaks of the 
(x, A) plot only depends on the energetic term. As we are 
considering a protein composed by identical domains, we 
expect the x(X) peaks to be identical if the temperature 
is sufficiently low. This pattern can be observed in Fig- 
ure ^2 where we compare two simulations performed at 
different temperatures, namely (3 = 2 and (3 — 8. 

D. Pulling rate effects 

Up to now, we discussed equilibrium stretching com- 
putations, i.e. we assumed an infinitely slow pulling. 
However, at the typical time scales of an AFM stretch- 
ing experiment, the polymer is pulled at a finite velocity 
and the unfolding is a non-equilibrium process, as tes- 
tified by the differences between the unfolding and the 
refolding force-extension profiles (see e.g. Figure IT^ . In 
fact, while the unfolding profile features the typical saw- 
tooth pattern, upon relaxation of the unfolded polymer, 
the trace exhibits no discontinuities that would indicate 
refolding. 

A consequence of the irreversibility of the stretching 
process is that the unfolding force depends on the pulling 
speed. Actually, if the loading rate rj — kfV (with kf be- 
ing the elastic constant and v the pulling velocity) is suf- 
ficiently small, thermal fluctuations are allowed enough 
time to overcome the energy barrier and the unfolding 
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FIG. 12: Monte Carlo computations of stretching and re- 
folding: the non-superposability of the profiles is a signature 
of the irreversibility of the process. Simulation parameters: 
N = 100; M = 3; L = 450; k x = 100; T = 1000; A = 1; 
9 = 0.3; K = 0.5; = 2. 



force will be low. 

Several experimental works 0, H3, El H2 reported 
a logarithmic dependence of the unfolding force on the 
loading rate in the case where a single energy barrier is 
present along the reaction path. The analytic expression 
of the relation between force and loading rate was de- 
rived [Tslls^lls^ l within the frame of Kramer's theory for 
a simple two-state model, by assuming that the external 
force reduces the height of the energy barrier, 



k B T ( KvAx \ 



(4) 



where ks is Boltzmann constant, T is the absolute tem- 
perature, Ax is the distance between the minimum cor- 
responding to the folded state and the activation barrier 
of the energy landscape, v is the pulling speed, K is the 
cantilever harmonic constant and ko is the spontaneous 
unfolding rate. 

Recently 19], it has been shown that the probability 
distribution of the force at various pulling rates does not 
follow the simple Bell law, requiring the full Kramer's 
theory and predicting small corrections to the logaritmic 
behavior of the most probable force. 

We limit here to a preliminary illustration of a series of 
Monte Carlo simulations with different number of steps 
T considering the pulling speed as being proportional to 
l/T. A more detailed analysis will be presented in a 
forthcoming paper. 

Figure FLU shows that, as the number of Monte-Carlo 
(MC) steps is increased, the (x — A) plot becomes closer 
and closer to the profile yielded by the equilibrium simu- 
lation. In particular, a small number of MC steps, i.e. a 
fast pulling, corresponds to high peaks of the (x — X) pro- 
file, whereas larger numbers of steps correspond to less 



FIG. 13: Comparison between the equilibrium computations 
and a set of Monte Carlo simulations with different numbers 
of MC steps. The peaks of the (x — A) profile become lower 
and lower and the MC simulations converge to the equilibrium 
scenario as the number of MC steps is increased. Simulation 
parameters: N = 10; M = 3; L = 70; k x = 100; A = 1; 
8 = 0.2; K = 0.1; = 2. Notice that, using a smaller value of 
N, the equilibrium profile is smoother than that in Figure |5] 
Figure rm and preceding ones. 
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FIG. 14: Statistical analysis of unfolding forces in 5 sets of 
100 independent Monte Carlo runs with different numbers of 
MC steps. Panel (a): histograms of spring extensions at the 
rupture point; panel (b): linear fit of the average rupture force 
as a function of the logarithm of the loading rate. Simulation 
parameters: N = 10; M = 3; L = 70; fc A = 100; A = 10; 
9 = 0.2; K = 0.1; = 2. 



and less pronounced peaks and thus, smaller unfolding 
forces. 

We performed a series of 100 independent MC runs for 
5 different values of T: 1 ■ 10 3 , 5 ■ 10 3 , 1 • 10 4 , 5 • 10 4 and 
1 • 10 5 . For each series of runs we built the histograms of 
the rupture spring elongations and we plotted the mean 
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value of the histogram as a function of the logarithm 
of the loading rate. Finally, the set of points thus ob- 
tained was linearly fitted. As shown in Figured! the his- 
tograms, shift to lower values as the number of MC steps 
increases. Figure 1141 also shows that the mean rupture 
forces computed from the histograms, feature an excel- 
lent linear correlation with the logarithm of the loading 
rate (correlation coefficient r c = 0.99). 

If Eq. 0] is explicitly rearranged so as to show the lin- 
ear dependence of (F max ) on log(Kv), the coefficients of 
the linear equation can be equated to the corresponding 
parameters 71 and 72 of the regression line (F max ) = 
71 log(Kv) +72, so as to build the following system of 
equations: 



k B T 
Ax 



Ax 
k k B T 



71 
72 



From the first equation of the system it is possible to 
compute the width of the activation barrier Ax = 11.17; 
this value can then be substituted into the second equa- 
tion so as to determine k = 2.93 • 10~ 7 . The spontaneous 
unfolding rate fco, on turn, is related to the height of the 
activation barrier for the unfolding process, 



fc 



: toe 



-AG'l/k B T 



where w, as explained by Kramer's theory, is the recipro- 
cal of a diffusive relaxation time. This simple computa- 
tion shows how stretching experiments and simulations 
provide easy access to important features of the free en- 
ergy landscape. In a forthcoming paper we aim at inves- 
tigating the relations relating the barrier width Ax and 
the spontaneous unfolding rate ko with molecular prop- 
erties such as the folding prize A, the threshold 9, the 
number N and the length M of the domains. 



E. High pulling rate limit 

Let us investigate the dependence of the Monte Carlo 
simulations on the number of Monte-Carlo steps T, that 
we can interpret as a measure of the pulling speed. In 
fact, in the limit of an extremely high pulling speed, we 
can assume that for each value of A the polymer can 
adopt just a single (or at most, a few) conformation, so 
that the entropic contribution can be neglected in the 
computations. 

The mechanism outlined in Section IIII Bl shows that 
the key event is the loss of the folding prize determining 
the complete extension of the protein domain and the 
subsequent relaxation of the spring. The problem thus 
arises to identify the factors affecting the probability of 
this crucial step. More specifically, suppose in a domain 
a number of contacts just below the threshold to keep the 
folding prize has been broken. What is the probability tt 
that one more contact will be cleaved within the next v 
steps ? 



In order to answer this question, it is necessary to com- 
pute the energy difference associated with the transition. 
Let x be the length of the spring when a number of con- 
tacts just below the threshold ON has been disrupted: 
the system retains the prize A. If one more contact is 
broken, the prize is lost and the spring will be corre- 
spondingly shortened. In particular, by assuming each 
contact breakdown to cause a unit increase in length of 
the polymer, the new length of the spring will be x — 1. 
The energy difference between the two states can thus be 
computed as 



\E = ( ^Kx 2 - Aj - \k(x - l) 2 = Kx 



A, 



where a term K/2 has been neglected. We can now com- 
pute the probability to destroy a contact in a single step 



1 



1 



P 



1 



-0AE 



I + e l3{A=Kx) ■ 



Finally, the probability to break a contact in v steps is 
given by the sum of a geometric progression 

V 

i=0 

The computations thus show that the domain unfolds 
with a probability depending on the prize, the cantilever 
stiffness, the temperature and the pulling velocity (that 
in our simulations is related to the number T of Monte 
Carlo steps). 
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FIG. 15: Effect of the simulation parameters on the height of 
the unfolding ramp. Unless otherwise indicated in the legend, 
the simulation parameters are: TV = 100; M = 1; L = 150; 
k x = 50; T = 100; A = 0.1; 9 = 0.2; K = 0.1; (3 = 2. 

In order to test our analytical treatment of the unfold- 
ing probability, we performed a series of MC simulations 
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differing from a reference one just for a parameter (Fig- 
ure !15|) . The key parameters characterizing the reference 
simulations are: T = 100, A = 0.1, K = 0.1, = 2. 
In agreement with the analytic computation, the simula- 
tions show that an increase in the number of Monte Carlo 
steps (T — 10 4 ), an increase in temperature (f3 — 1.5), 
a decrease of the folding prize (A = 8.5 x 10~ 2 ), and an 
increase of the spring constant (K — 1.25 x 10 _1 ), all 
result in an increase of the unfolding probability, thus 
reducing the height of the peak in the (x — A)-plot. 

IV. CONCLUSIONS 

We developed an extremely simplified model of poly- 
mer stretching in which the molecule is portrayed as a 
series of modules, represented as an array of contacts, and 
a harmonic spring (the cantilever). The chemical in- 
teractions stabilizing the native conformation are sim- 
ply modeled as a folding prize gained by domains where 
the fraction of folded monomers is above a pre-chosen 
threshold. Our model is consistent with recent findings 
in Refs [2lL I22I I23I ] , showing that the unraveling of the 
titin Immunoglobulin domain occurs very rapidly only 
after the breakdown of a critical number of key hydrogen- 
bonds. The attribution of a folding prize when a thresh- 
old value of the fraction of native contacts is crossed, also 
makes our approach equivalent to a Go-model with force 
rescaling |26j. However, our model is thus significantly 
simpler than other models commonly used to study me- 
chanical unfolding such as the WLC and FJC models 
and the detailed all-atom, topological and united-residue 
models employed for steered molecular dynamics simu- 
lations. Yet, our model is detailed enough to reproduce 
many qualitative features of the force-extension profiles 
recorded in AFM experiments. In particular, our model 
correctly reproduces the typical saw-tooth pattern with 
peaks characterized by a height dependent on the folding 
prize, the temperature and the pulling velocity. 

In our study, a particular attention was paid to the 
relation between the heights of successive peaks in the 
force-extension plot. This point is quite intriguing as ex- 
perimental data show that the force of unfolding tends 
to increase with each unfolding event, suggesting that 
the protein domains unfold following an increasing order 
of mechanical stability j^Ej. A possible explanation sug- 
gested by the literature is that the domains of the engi- 
neered modular constructs used in AFM experiments are 
not identical but just structurally similar 15] . This ar- 
gument may be correct in the case of real proteins where 
the differences in mechanical stability of the tandem re- 
peats of the construct may arise from their different po- 
sition in the tridimensional structure of the molecule. In 



a computer simulation, however, the domains are all per- 
fectly identical and the explanation for the staircase pat- 
tern of unfolding peaks must be sought elsewhere. This 
behavior has been ascribed to an independent breaking 
probability of each domain. This hypothesis is however 
not consistent with the fact that the cantilever couples 
the fluctuations of all domains. Our unified treatment 
allows to keep into consideration all contributions, and 
to obtain the correct equilibrium curves, that generally 
do not exhibit this effect. However, a detailed study of 
out-of-cquilibrium extension curves, shows that this be- 
havior is consistent with a finite pulling speed, i.e. it can 
be ascribed to a dynamical source. 

Our results appear to be in agreement with recent find- 
ings by Cieplak 36] in steered molecular dynamics simu- 
lations of titin and calmodulin unfolding using a Go-like 
model. In particular, it was found that an increase in ther- 
mal fluctuations results in a lowering of the force peaks 
and in their earlier occurrence during stretching. 

An interesting feature of AFM stretching experiments 
is that the mean unfolding force is a logarithmic function 
of the pulling rate. This pattern reflects the role of ther- 
mal fluctuations: if the pulling speed is sufficiently low, 
then there will be enough time for thermal fluctuations 
to drive the polymer over the free energy barrier and the 
unfolding force will be accordingly low. The ability to 
reproduce this pattern constitutes a benchmark for any 
stretching model. In order to test our model, we per- 
formed a series of Monte Carlo simulations with different 
numbers T of MC-steps and regarding the pulling veloc- 
ity to be proportional to 1/T. Our toy-model, despite its 
extreme simplicity, correctly reproduced the linear de- 
pendence of the mean unfolding force on the logarithm 
of the pulling rate. 

In summary, our simple model including harmonic and 
long-range energy contributions, is not only capable of re- 
producing the force-extension saw-tooth pattern, but it 
also yields force spectra displaying the correct logarith- 
mic dependence of (F max ) on Kv reported in experimen- 
tal works. Finally, our model provides a simple expla- 
nation of the influence of each unfolding event on the 
following ones, based on the role of fluctuations. Our 
work thus shows how minimal models can be valuable 
tools even in the study of complex molecular systems. 
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